This is code used to make plots based on output from the “IndividualPondModels_Forecasts85” for boththe CR and YF Read in the output files for each pond/RCP combination then work on delta figures and figures like Hamlet et al. 2020
rm(list = ls())
library(tidyverse)
library(dplyr)
library(rjags)
library(ggplot2)
set.seed(1)
Starting with the Yakutat Forelands
RCP 85
Mod_YF85 <- read.csv("Mod_YF85.csv", header = T)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_YF85$pond_name)
# Initialize an empty dataframe to store results
avgMed_2020 <- data.frame(Pond = pond_names, AvgMed_2020 = numeric(length(pond_names)))
# Loop through each pond name
for (i in seq_along(pond_names)) {
pond_name <- pond_names[i]
# Filter data for the current pond and for dates before January 1, 2021
filtered_data <- Mod_YF85 %>%
filter(pond_name == !!pond_name, as.Date(time) < as.Date("2021-01-01"))
# Calculate the mean of 'med' (which corresponds to X50%)
avg_med <- mean(filtered_data$med, na.rm = TRUE)
# Store the result in the dataframe
avgMed_2020[i, "AvgMed_2020"] <- avg_med
}
# Print the dataframe
print(avgMed_2020)
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_YF85$pond_name)
# Initialize a dataframe to store results with years as columns
avgMed_2030 <- data.frame(Pond = pond_names)
# Loop through each year from 2030 to 2039
for (yr in 2030:2039) {
# Calculate average medians for each year and store as a new column
avgMed_col <- sapply(pond_names, function(pond_name) {
yearly_data <- Mod_YF85 %>%
filter(pond_name == !!pond_name, year(as.Date(time)) == yr)
avg_med <- mean(yearly_data$med, na.rm = TRUE)
return(avg_med)
})
# Add the results as a new column to the dataframe
avgMed_2030[[as.character(yr)]] <- avgMed_col
}
# Print the final dataframe
print(avgMed_2030)
# Convert to long format for ggplot2
avgMed_2030_long <- avgMed_2030 %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed")
# Plot the data
ggplot(avgMed_2030_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond from 2030 to 2039",
x = "Year",
y = "Average Median") +
theme_minimal()
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_YF85$pond_name)
# Initialize a dataframe to store results with years as columns
avgMed_2060 <- data.frame(Pond = pond_names)
# Loop through each year from 2060 to 2069
for (yr in 2060:2069) {
# Calculate average medians for each year and store as a new column
avgMed_col <- sapply(pond_names, function(pond_name) {
yearly_data <- Mod_YF85 %>%
filter(pond_name == !!pond_name, year(as.Date(time)) == yr)
avg_med <- mean(yearly_data$med, na.rm = TRUE)
return(avg_med)
})
# Add the results as a new column to the dataframe
avgMed_2060[[as.character(yr)]] <- avgMed_col
}
# Print the final dataframe
print(avgMed_2060)
# Convert to long format for ggplot2
avgMed_2060_long <- avgMed_2060 %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed")
# Plot the data
ggplot(avgMed_2060_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond from 2060 to 2069",
x = "Year",
y = "Average Median") +
theme_minimal()
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_YF85$pond_name)
# Initialize a dataframe to store results with years as columns
avgMed_2090 <- data.frame(Pond = pond_names)
# Loop through each year from 2090 to 2099
for (yr in 2090:2099) {
# Calculate average medians for each year and store as a new column
avgMed_col <- sapply(pond_names, function(pond_name) {
yearly_data <- Mod_YF85 %>%
filter(pond_name == !!pond_name, year(as.Date(time)) == yr)
avg_med <- mean(yearly_data$med, na.rm = TRUE)
return(avg_med)
})
# Add the results as a new column to the dataframe
avgMed_2090[[as.character(yr)]] <- avgMed_col
}
# Print the final dataframe
print(avgMed_2090)
# Convert to long format for ggplot2
avgMed_2090_long <- avgMed_2090 %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed")
# Plot the data
ggplot(avgMed_2090_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond from 2090 to 2099",
x = "Year",
y = "Average Median") +
theme_minimal()
# Rename columns to reflect only the years
avgMed_2020 <- avgMed_2020 %>%
rename(AvgMed_2020 = AvgMed_2020)
avgMed_2030 <- avgMed_2030 %>%
rename(`2030` = `2030`)
avgMed_2060 <- avgMed_2060 %>%
rename(`2060` = `2060`)
avgMed_2090 <- avgMed_2090 %>%
rename(`2090` = `2090`)
# Combine the data frames
combined_avgMed <- avgMed_2020 %>%
left_join(avgMed_2030, by = "Pond") %>%
left_join(avgMed_2060, by = "Pond") %>%
left_join(avgMed_2090, by = "Pond")
# Print the final dataframe
print(combined_avgMed)
# Convert to long format for ggplot2, excluding AvgMed_2020
avgMed_long <- combined_avgMed %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed") %>%
filter(Year %in% as.character(2021:2199))
# Plot the data
ggplot(avgMed_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond in the YF - RCP 8.5",
x = "Year",
y = "Average Median") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),# Rotate x-axis labels by 45 degrees
plot.title = element_text(hjust = 0.5))
library(dplyr)
# Assuming combined_avgMed is loaded and structure verified
Delta_combined_avgMed <- combined_avgMed # Use your actual dataframe
# Columns to calculate differences for
years <- c("2030", "2031", "2032", "2033", "2034",
"2035", "2036", "2037", "2038", "2039",
"2060", "2061", "2062", "2063", "2064",
"2065", "2066", "2067", "2068", "2069",
"2090", "2091", "2092", "2093", "2094",
"2095", "2096", "2097", "2098", "2099")
# Loop through each year and calculate the difference
for (year in years) {
# Create a new column for the difference
diff_col_name <- paste0("diff_", year)
Delta_combined_avgMed[[diff_col_name]] <- Delta_combined_avgMed[[year]] - Delta_combined_avgMed$`AvgMed_2020`
}
# Display the structure of Delta_combined_avgMed to verify
str(Delta_combined_avgMed)
'data.frame': 11 obs. of 62 variables:
$ Pond : chr "MP1" "UBP3" "UBP4" "MP3" ...
$ AvgMed_2020: num 7.45 7.13 7.87 5.47 6.77 ...
$ 2030 : num 9.42 9.09 9.85 7.41 8.74 ...
$ 2031 : num 6.82 6.5 7.27 4.82 6.14 ...
$ 2032 : num 8.46 8.13 8.88 6.44 7.78 ...
$ 2033 : num 9.43 9.1 9.85 7.42 8.72 ...
$ 2034 : num 7.97 7.63 8.38 5.95 7.26 ...
$ 2035 : num 8.89 8.55 9.3 6.87 8.19 ...
$ 2036 : num 8.66 8.33 9.09 6.65 7.96 ...
$ 2037 : num 10.12 9.8 10.54 8.11 9.43 ...
$ 2038 : num 9.55 9.24 9.98 7.54 8.87 ...
$ 2039 : num 8.12 7.81 8.55 6.11 7.44 ...
$ 2060 : num 10.86 10.53 11.28 8.84 10.17 ...
$ 2061 : num 11.09 10.78 11.53 9.09 10.4 ...
$ 2062 : num 12.5 12.2 13 10.6 11.9 ...
$ 2063 : num 14.2 13.9 14.7 12.2 13.5 ...
$ 2064 : num 11.6 11.28 12.04 9.59 10.92 ...
$ 2065 : num 10.98 10.66 11.41 8.97 10.29 ...
$ 2066 : num 12.3 11.9 12.7 10.2 11.6 ...
$ 2067 : num 11.16 10.82 11.58 9.15 10.47 ...
$ 2068 : num 11.29 10.98 11.74 9.29 10.6 ...
$ 2069 : num 12.3 12 12.7 10.3 11.6 ...
$ 2090 : num 14.2 13.9 14.6 12.2 13.5 ...
$ 2091 : num 14.8 14.5 15.2 12.8 14.1 ...
$ 2092 : num 14.3 13.9 14.7 12.3 13.6 ...
$ 2093 : num 14.5 14.2 14.9 12.5 13.8 ...
$ 2094 : num 13.9 13.6 14.4 11.9 13.2 ...
$ 2095 : num 14.4 14.1 14.8 12.3 13.7 ...
$ 2096 : num 15.2 14.9 15.6 13.2 14.5 ...
$ 2097 : num 14.7 14.3 15.1 12.7 14 ...
$ 2098 : num 14.9 14.6 15.3 12.9 14.2 ...
$ 2099 : num 15.3 15 15.7 13.3 14.6 ...
$ diff_2030 : num 1.97 1.96 1.98 1.95 1.97 ...
$ diff_2031 : num -0.625 -0.633 -0.607 -0.647 -0.629 ...
$ diff_2032 : num 1.01 1.006 1.009 0.972 1.015 ...
$ diff_2033 : num 1.98 1.98 1.98 1.95 1.96 ...
$ diff_2034 : num 0.516 0.502 0.509 0.484 0.493 ...
$ diff_2035 : num 1.44 1.43 1.43 1.4 1.42 ...
$ diff_2036 : num 1.21 1.2 1.21 1.18 1.2 ...
$ diff_2037 : num 2.67 2.67 2.67 2.64 2.66 ...
$ diff_2038 : num 2.1 2.11 2.1 2.07 2.1 ...
$ diff_2039 : num 0.675 0.679 0.679 0.639 0.671 ...
$ diff_2060 : num 3.41 3.41 3.41 3.38 3.4 ...
$ diff_2061 : num 3.64 3.65 3.66 3.62 3.63 ...
$ diff_2062 : num 5.1 5.09 5.11 5.1 5.1 ...
$ diff_2063 : num 6.78 6.78 6.78 6.76 6.76 ...
$ diff_2064 : num 4.15 4.15 4.17 4.12 4.15 ...
$ diff_2065 : num 3.53 3.54 3.54 3.51 3.52 ...
$ diff_2066 : num 4.81 4.82 4.82 4.78 4.82 ...
$ diff_2067 : num 3.71 3.7 3.71 3.68 3.7 ...
$ diff_2068 : num 3.84 3.85 3.86 3.82 3.83 ...
$ diff_2069 : num 4.84 4.82 4.84 4.81 4.83 ...
$ diff_2090 : num 6.75 6.75 6.77 6.71 6.74 ...
$ diff_2091 : num 7.39 7.38 7.37 7.35 7.35 ...
$ diff_2092 : num 6.82 6.82 6.82 6.79 6.82 ...
$ diff_2093 : num 7.04 7.02 7.03 7.01 7.04 ...
$ diff_2094 : num 6.49 6.49 6.5 6.47 6.48 ...
$ diff_2095 : num 6.92 6.93 6.91 6.88 6.9 ...
$ diff_2096 : num 7.75 7.76 7.76 7.72 7.75 ...
$ diff_2097 : num 7.22 7.22 7.24 7.2 7.22 ...
$ diff_2098 : num 7.45 7.44 7.46 7.42 7.44 ...
$ diff_2099 : num 7.85 7.84 7.85 7.83 7.84 ...
library(matrixStats)
# Create a new dataframe with the calculated means and standard deviations
Delta_YF85 <- Delta_combined_avgMed %>%
mutate(
Mean_2030_2039 = rowMeans(select(., starts_with("diff_203")), na.rm = TRUE), # Calculate row-wise mean for columns starting with "diff_203"
SD_2030_2039 = rowSds(as.matrix(select(., starts_with("diff_203"))), na.rm = TRUE), # Calculate row-wise SD for columns starting with "diff_203"
Mean_2060_2069 = rowMeans(select(., starts_with("diff_206")), na.rm = TRUE), # Calculate row-wise mean for columns starting with "diff_206"
SD_2060_2069 = rowSds(as.matrix(select(., starts_with("diff_206"))), na.rm = TRUE), # Calculate row-wise SD for columns starting with "diff_206"
Mean_2090_2099 = rowMeans(select(., starts_with("diff_209")), na.rm = TRUE), # Calculate row-wise mean for columns starting with "diff_209"
SD_2090_2099 = rowSds(as.matrix(select(., starts_with("diff_209"))), na.rm = TRUE) # Calculate row-wise SD for columns starting with "diff_209"
)
# Print the structure of the new dataframe
str(Delta_YF85)
'data.frame': 11 obs. of 68 variables:
$ Pond : chr "MP1" "UBP3" "UBP4" "MP3" ...
$ AvgMed_2020 : num 7.45 7.13 7.87 5.47 6.77 ...
$ 2030 : num 9.42 9.09 9.85 7.41 8.74 ...
$ 2031 : num 6.82 6.5 7.27 4.82 6.14 ...
$ 2032 : num 8.46 8.13 8.88 6.44 7.78 ...
$ 2033 : num 9.43 9.1 9.85 7.42 8.72 ...
$ 2034 : num 7.97 7.63 8.38 5.95 7.26 ...
$ 2035 : num 8.89 8.55 9.3 6.87 8.19 ...
$ 2036 : num 8.66 8.33 9.09 6.65 7.96 ...
$ 2037 : num 10.12 9.8 10.54 8.11 9.43 ...
$ 2038 : num 9.55 9.24 9.98 7.54 8.87 ...
$ 2039 : num 8.12 7.81 8.55 6.11 7.44 ...
$ 2060 : num 10.86 10.53 11.28 8.84 10.17 ...
$ 2061 : num 11.09 10.78 11.53 9.09 10.4 ...
$ 2062 : num 12.5 12.2 13 10.6 11.9 ...
$ 2063 : num 14.2 13.9 14.7 12.2 13.5 ...
$ 2064 : num 11.6 11.28 12.04 9.59 10.92 ...
$ 2065 : num 10.98 10.66 11.41 8.97 10.29 ...
$ 2066 : num 12.3 11.9 12.7 10.2 11.6 ...
$ 2067 : num 11.16 10.82 11.58 9.15 10.47 ...
$ 2068 : num 11.29 10.98 11.74 9.29 10.6 ...
$ 2069 : num 12.3 12 12.7 10.3 11.6 ...
$ 2090 : num 14.2 13.9 14.6 12.2 13.5 ...
$ 2091 : num 14.8 14.5 15.2 12.8 14.1 ...
$ 2092 : num 14.3 13.9 14.7 12.3 13.6 ...
$ 2093 : num 14.5 14.2 14.9 12.5 13.8 ...
$ 2094 : num 13.9 13.6 14.4 11.9 13.2 ...
$ 2095 : num 14.4 14.1 14.8 12.3 13.7 ...
$ 2096 : num 15.2 14.9 15.6 13.2 14.5 ...
$ 2097 : num 14.7 14.3 15.1 12.7 14 ...
$ 2098 : num 14.9 14.6 15.3 12.9 14.2 ...
$ 2099 : num 15.3 15 15.7 13.3 14.6 ...
$ diff_2030 : num 1.97 1.96 1.98 1.95 1.97 ...
$ diff_2031 : num -0.625 -0.633 -0.607 -0.647 -0.629 ...
$ diff_2032 : num 1.01 1.006 1.009 0.972 1.015 ...
$ diff_2033 : num 1.98 1.98 1.98 1.95 1.96 ...
$ diff_2034 : num 0.516 0.502 0.509 0.484 0.493 ...
$ diff_2035 : num 1.44 1.43 1.43 1.4 1.42 ...
$ diff_2036 : num 1.21 1.2 1.21 1.18 1.2 ...
$ diff_2037 : num 2.67 2.67 2.67 2.64 2.66 ...
$ diff_2038 : num 2.1 2.11 2.1 2.07 2.1 ...
$ diff_2039 : num 0.675 0.679 0.679 0.639 0.671 ...
$ diff_2060 : num 3.41 3.41 3.41 3.38 3.4 ...
$ diff_2061 : num 3.64 3.65 3.66 3.62 3.63 ...
$ diff_2062 : num 5.1 5.09 5.11 5.1 5.1 ...
$ diff_2063 : num 6.78 6.78 6.78 6.76 6.76 ...
$ diff_2064 : num 4.15 4.15 4.17 4.12 4.15 ...
$ diff_2065 : num 3.53 3.54 3.54 3.51 3.52 ...
$ diff_2066 : num 4.81 4.82 4.82 4.78 4.82 ...
$ diff_2067 : num 3.71 3.7 3.71 3.68 3.7 ...
$ diff_2068 : num 3.84 3.85 3.86 3.82 3.83 ...
$ diff_2069 : num 4.84 4.82 4.84 4.81 4.83 ...
$ diff_2090 : num 6.75 6.75 6.77 6.71 6.74 ...
$ diff_2091 : num 7.39 7.38 7.37 7.35 7.35 ...
$ diff_2092 : num 6.82 6.82 6.82 6.79 6.82 ...
$ diff_2093 : num 7.04 7.02 7.03 7.01 7.04 ...
$ diff_2094 : num 6.49 6.49 6.5 6.47 6.48 ...
$ diff_2095 : num 6.92 6.93 6.91 6.88 6.9 ...
$ diff_2096 : num 7.75 7.76 7.76 7.72 7.75 ...
$ diff_2097 : num 7.22 7.22 7.24 7.2 7.22 ...
$ diff_2098 : num 7.45 7.44 7.46 7.42 7.44 ...
$ diff_2099 : num 7.85 7.84 7.85 7.83 7.84 ...
$ Mean_2030_2039: num 1.29 1.29 1.3 1.26 1.29 ...
$ SD_2030_2039 : num 0.957 0.96 0.955 0.957 0.958 ...
$ Mean_2060_2069: num 4.38 4.38 4.39 4.36 4.38 ...
$ SD_2060_2069 : num 1.04 1.03 1.03 1.04 1.03 ...
$ Mean_2090_2099: num 7.17 7.17 7.17 7.14 7.16 ...
$ SD_2090_2099 : num 0.443 0.44 0.441 0.443 0.444 ...
# Optionally, view the first few rows to check the new columns
head(Delta_YF85)
# Load necessary libraries
library(dplyr)
# Select and rename columns from Delta_YF85
YF85_decadal <- Delta_YF85 %>%
select(Pond, starts_with("Mean_"), starts_with("SD_"))
# Rename columns without the 'Y' prefix for means
names(YF85_decadal) <- sub("^Mean_", "Mean_", names(YF85_decadal))
# Rename columns without the 'Y' prefix for standard deviations
names(YF85_decadal) <- sub("^SD_", "SD_", names(YF85_decadal))
# Print the structure of the updated dataframe
str(YF85_decadal)
'data.frame': 11 obs. of 7 variables:
$ Pond : chr "MP1" "UBP3" "UBP4" "MP3" ...
$ Mean_2030_2039: num 1.29 1.29 1.3 1.26 1.29 ...
$ Mean_2060_2069: num 4.38 4.38 4.39 4.36 4.38 ...
$ Mean_2090_2099: num 7.17 7.17 7.17 7.14 7.16 ...
$ SD_2030_2039 : num 0.957 0.96 0.955 0.957 0.958 ...
$ SD_2060_2069 : num 1.04 1.03 1.03 1.04 1.03 ...
$ SD_2090_2099 : num 0.443 0.44 0.441 0.443 0.444 ...
# Optionally, view the first few rows to check the new columns
head(YF85_decadal)
# Load required packages
library(ggplot2)
library(gridExtra)
library(ggpubr)
# Create a grouped bar plot with error bars for Mean 2030
plot_2030 <- ggplot(YF85_decadal, aes(x = Pond, y = Mean_2030_2039)) +
geom_bar(stat = "identity", fill = "#FFD92F", color = "black") +
geom_errorbar(aes(ymin = Mean_2030_2039 - SD_2030_2039, ymax = Mean_2030_2039 + SD_2030_2039),
width = 0.4, # Adjust the width of error bars as needed
position = position_dodge(width = 0.9)) +
labs(x = "", y = "Delta Temp (C)") +
ylim(-1, 5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) # Centered title
plot_2030
# Create a grouped bar plot with error bars for Mean 2060
plot_2060 <- ggplot(YF85_decadal, aes(x = Pond, y = Mean_2060_2069)) +
geom_bar(stat = "identity", fill = "#F46D43", color = "black") +
geom_errorbar(aes(ymin = Mean_2060_2069 - SD_2060_2069, ymax = Mean_2060_2069 + SD_2060_2069),
width = 0.4, # Adjust the width of error bars as needed
position = position_dodge(width = 0.9)) +
labs(x = "", y = "") + # No y-axis label
ylim(-1, 5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) + # Centered title
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank()) # Remove Y-axis ticks
plot_2060
# Create a grouped bar plot with error bars for Mean 2090
plot_2090 <- ggplot(YF85_decadal, aes(x = Pond, y = Mean_2090_2099)) +
geom_bar(stat = "identity", fill = "#D73027", color = "black") +
geom_errorbar(aes(ymin = Mean_2090_2099 - SD_2090_2099, ymax = Mean_2090_2099 + SD_2090_2099),
width = 0.4, # Adjust the width of error bars as needed
position = position_dodge(width = 0.9)) +
labs(x = "", y = "") + # No y-axis label
ylim(-1, 5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) + # Centered title
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank()) # Remove Y-axis ticks
plot_2090
Warning: Removed 11 rows containing missing values or values outside the scale range (`geom_bar()`).
# Combine plots into a panel plot
Delta_YF85_plot <- ggarrange(plot_2030, plot_2060, plot_2090, ncol = 3)
Warning: Removed 11 rows containing missing values or values outside the scale range (`geom_bar()`).
# Save the combined plot
ggsave("Corr_DeltaYF85.png", plot = Delta_YF85_plot, width = 8, height = 3)
Delta_YF85
# Load required libraries
library(dplyr)
library(broom)
library(tidyr)
# Step 1: Prepare the data
data_long <- Delta_YF85 %>%
select(Pond, matches("^diff_203[0-9]$")) %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "Difference")
boxplot(data_long$Difference ~ data_long$Pond)
# Step 2: Conduct ANOVA by grouping by Pond
res_YF85 <- aov(Difference ~ Pond, data = data_long)
summary(res_YF85)
Df Sum Sq Mean Sq F value Pr(>F)
Pond 10 0.01 0.0007 0.001 1
Residuals 99 90.74 0.9165
# Step 3: Tukey HSD
post_test_YF85 <- TukeyHSD(res_YF85)
post_test_YF85
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Difference ~ Pond, data = data_long)
$Pond
diff lwr upr p adj
MP3-MP1 -0.0309536340 -1.441778 1.379871 1
MP5-MP1 -0.0084655775 -1.419290 1.402359 1
MP8-MP1 -0.0070506871 -1.417875 1.403774 1
PL1-MP1 -0.0053581246 -1.416183 1.405467 1
PL2-MP1 -0.0077879959 -1.418613 1.403037 1
PL3-MP1 -0.0104036062 -1.421228 1.400421 1
UBP1-MP1 -0.0069344459 -1.417759 1.403890 1
UBP2-MP1 -0.0086695815 -1.419494 1.402155 1
UBP3-MP1 -0.0046855885 -1.415510 1.406139 1
UBP4-MP1 0.0018912822 -1.408933 1.412716 1
MP5-MP3 0.0224880565 -1.388337 1.433313 1
MP8-MP3 0.0239029468 -1.386922 1.434728 1
PL1-MP3 0.0255955093 -1.385229 1.436420 1
PL2-MP3 0.0231656381 -1.387659 1.433990 1
PL3-MP3 0.0205500278 -1.390275 1.431375 1
UBP1-MP3 0.0240191880 -1.386805 1.434844 1
UBP2-MP3 0.0222840525 -1.388541 1.433109 1
UBP3-MP3 0.0262680454 -1.384557 1.437093 1
UBP4-MP3 0.0328449161 -1.377980 1.443670 1
MP8-MP5 0.0014148904 -1.409410 1.412240 1
PL1-MP5 0.0031074529 -1.407717 1.413932 1
PL2-MP5 0.0006775816 -1.410147 1.411502 1
PL3-MP5 -0.0019380287 -1.412763 1.408887 1
UBP1-MP5 0.0015311316 -1.409294 1.412356 1
UBP2-MP5 -0.0002040040 -1.411029 1.410621 1
UBP3-MP5 0.0037799890 -1.407045 1.414605 1
UBP4-MP5 0.0103568597 -1.400468 1.421181 1
PL1-MP8 0.0016925625 -1.409132 1.412517 1
PL2-MP8 -0.0007373087 -1.411562 1.410087 1
PL3-MP8 -0.0033529190 -1.414178 1.407472 1
UBP1-MP8 0.0001162412 -1.410708 1.410941 1
UBP2-MP8 -0.0016188943 -1.412444 1.409206 1
UBP3-MP8 0.0023650986 -1.408460 1.413190 1
UBP4-MP8 0.0089419693 -1.401883 1.419767 1
PL2-PL1 -0.0024298713 -1.413255 1.408395 1
PL3-PL1 -0.0050454816 -1.415870 1.405779 1
UBP1-PL1 -0.0015763213 -1.412401 1.409248 1
UBP2-PL1 -0.0033114569 -1.414136 1.407513 1
UBP3-PL1 0.0006725361 -1.410152 1.411497 1
UBP4-PL1 0.0072494068 -1.403575 1.418074 1
PL3-PL2 -0.0026156103 -1.413440 1.408209 1
UBP1-PL2 0.0008535499 -1.409971 1.411678 1
UBP2-PL2 -0.0008815856 -1.411706 1.409943 1
UBP3-PL2 0.0031024074 -1.407722 1.413927 1
UBP4-PL2 0.0096792781 -1.401145 1.420504 1
UBP1-PL3 0.0034691603 -1.407355 1.414294 1
UBP2-PL3 0.0017340247 -1.409091 1.412559 1
UBP3-PL3 0.0057180177 -1.405107 1.416543 1
UBP4-PL3 0.0122948884 -1.398530 1.423120 1
UBP2-UBP1 -0.0017351355 -1.412560 1.409090 1
UBP3-UBP1 0.0022488574 -1.408576 1.413073 1
UBP4-UBP1 0.0088257281 -1.401999 1.419650 1
UBP3-UBP2 0.0039839930 -1.406841 1.414809 1
UBP4-UBP2 0.0105608637 -1.400264 1.421385 1
UBP4-UBP3 0.0065768707 -1.404248 1.417402 1
Delta_YF85
# Load required libraries
library(dplyr)
library(broom)
# Step 1: Prepare the data
data_long <- Delta_YF85 %>%
select(Pond, matches("^diff_206[0-9]$")) %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "Difference")
boxplot(data_long$Difference ~ data_long$Pond)
# Step 2: Conduct ANOVA by grouping by Pond
res_YF85 <- aov(Difference ~ Pond, data = data_long)
summary(res_YF85)
Df Sum Sq Mean Sq F value Pr(>F)
Pond 10 0.01 0.0006 0.001 1
Residuals 99 106.35 1.0742
# Step 3: Tukey HSD
post_test_YF85 <- TukeyHSD(res_YF85)
post_test_YF85
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Difference ~ Pond, data = data_long)
$Pond
diff lwr upr p adj
MP3-MP1 -0.0239143576 -1.551308 1.503479 1
MP5-MP1 -0.0062309085 -1.533624 1.521162 1
MP8-MP1 -0.0006772301 -1.528070 1.526716 1
PL1-MP1 -0.0037509302 -1.531144 1.523642 1
PL2-MP1 -0.0091502094 -1.536543 1.518243 1
PL3-MP1 -0.0023098736 -1.529703 1.525083 1
UBP1-MP1 -0.0077866481 -1.535180 1.519607 1
UBP2-MP1 -0.0063965106 -1.533790 1.520997 1
UBP3-MP1 -0.0014816268 -1.528875 1.525912 1
UBP4-MP1 0.0080534452 -1.519340 1.535447 1
MP5-MP3 0.0176834491 -1.509710 1.545077 1
MP8-MP3 0.0232371275 -1.504156 1.550630 1
PL1-MP3 0.0201634274 -1.507230 1.547557 1
PL2-MP3 0.0147641481 -1.512629 1.542157 1
PL3-MP3 0.0216044840 -1.505789 1.548998 1
UBP1-MP3 0.0161277094 -1.511266 1.543521 1
UBP2-MP3 0.0175178469 -1.509875 1.544911 1
UBP3-MP3 0.0224327308 -1.504961 1.549826 1
UBP4-MP3 0.0319678027 -1.495425 1.559361 1
MP8-MP5 0.0055536784 -1.521840 1.532947 1
PL1-MP5 0.0024799783 -1.524913 1.529873 1
PL2-MP5 -0.0029193010 -1.530313 1.524474 1
PL3-MP5 0.0039210349 -1.523472 1.531314 1
UBP1-MP5 -0.0015557397 -1.528949 1.525838 1
UBP2-MP5 -0.0001656022 -1.527559 1.527228 1
UBP3-MP5 0.0047492817 -1.522644 1.532143 1
UBP4-MP5 0.0142843536 -1.513109 1.541678 1
PL1-MP8 -0.0030737001 -1.530467 1.524320 1
PL2-MP8 -0.0084729794 -1.535866 1.518920 1
PL3-MP8 -0.0016326435 -1.529026 1.525761 1
UBP1-MP8 -0.0071094181 -1.534503 1.520284 1
UBP2-MP8 -0.0057192806 -1.533113 1.521674 1
UBP3-MP8 -0.0008043967 -1.528198 1.526589 1
UBP4-MP8 0.0087306752 -1.518663 1.536124 1
PL2-PL1 -0.0053992793 -1.532793 1.521994 1
PL3-PL1 0.0014410566 -1.525952 1.528834 1
UBP1-PL1 -0.0040357180 -1.531429 1.523358 1
UBP2-PL1 -0.0026455805 -1.530039 1.524748 1
UBP3-PL1 0.0022693034 -1.525124 1.529663 1
UBP4-PL1 0.0118043753 -1.515589 1.539198 1
PL3-PL2 0.0068403358 -1.520553 1.534234 1
UBP1-PL2 0.0013635613 -1.526030 1.528757 1
UBP2-PL2 0.0027536988 -1.524640 1.530147 1
UBP3-PL2 0.0076685827 -1.519725 1.535062 1
UBP4-PL2 0.0172036546 -1.510190 1.544597 1
UBP1-PL3 -0.0054767745 -1.532870 1.521916 1
UBP2-PL3 -0.0040866370 -1.531480 1.523307 1
UBP3-PL3 0.0008282468 -1.526565 1.528221 1
UBP4-PL3 0.0103633187 -1.517030 1.537757 1
UBP2-UBP1 0.0013901375 -1.526003 1.528783 1
UBP3-UBP1 0.0063050214 -1.521088 1.533698 1
UBP4-UBP1 0.0158400933 -1.511553 1.543233 1
UBP3-UBP2 0.0049148839 -1.522478 1.532308 1
UBP4-UBP2 0.0144499558 -1.512943 1.541843 1
UBP4-UBP3 0.0095350719 -1.517858 1.536928 1
Delta_YF85
# Load required libraries
library(dplyr)
library(broom)
# Step 1: Prepare the data
data_long <- Delta_YF85 %>%
select(Pond, matches("^diff_209[0-9]$")) %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "Difference")
boxplot(data_long$Difference ~ data_long$Pond)
# Step 2: Conduct ANOVA by grouping by Pond
res_YF85 <- aov(Difference ~ Pond, data = data_long)
summary(res_YF85)
Df Sum Sq Mean Sq F value Pr(>F)
Pond 10 0.008 0.00081 0.004 1
Residuals 99 19.461 0.19657
# Step 3: Tukey HSD
post_test_YF85 <- TukeyHSD(res_YF85)
post_test_YF85
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Difference ~ Pond, data = data_long)
$Pond
diff lwr upr p adj
MP3-MP1 -3.113857e-02 -0.6845202 0.6222431 1
MP5-MP1 -9.882966e-03 -0.6632646 0.6434987 1
MP8-MP1 -1.672228e-03 -0.6550539 0.6517094 1
PL1-MP1 -4.243460e-03 -0.6576251 0.6491382 1
PL2-MP1 -6.484398e-03 -0.6598660 0.6468972 1
PL3-MP1 -9.998092e-03 -0.6633797 0.6433835 1
UBP1-MP1 -1.006286e-02 -0.6634445 0.6433188 1
UBP2-MP1 -6.892669e-03 -0.6602743 0.6464890 1
UBP3-MP1 -2.228373e-03 -0.6556100 0.6511533 1
UBP4-MP1 2.448355e-03 -0.6509333 0.6558300 1
MP5-MP3 2.125560e-02 -0.6321260 0.6746372 1
MP8-MP3 2.946634e-02 -0.6239153 0.6828480 1
PL1-MP3 2.689511e-02 -0.6264865 0.6802767 1
PL2-MP3 2.465417e-02 -0.6287275 0.6780358 1
PL3-MP3 2.114047e-02 -0.6322412 0.6745221 1
UBP1-MP3 2.107571e-02 -0.6323059 0.6744573 1
UBP2-MP3 2.424590e-02 -0.6291357 0.6776275 1
UBP3-MP3 2.891019e-02 -0.6244714 0.6822918 1
UBP4-MP3 3.358692e-02 -0.6197947 0.6869686 1
MP8-MP5 8.210738e-03 -0.6451709 0.6615924 1
PL1-MP5 5.639506e-03 -0.6477421 0.6590211 1
PL2-MP5 3.398567e-03 -0.6499831 0.6567802 1
PL3-MP5 -1.151261e-04 -0.6534968 0.6532665 1
UBP1-MP5 -1.798952e-04 -0.6535615 0.6532017 1
UBP2-MP5 2.990296e-03 -0.6503913 0.6563719 1
UBP3-MP5 7.654593e-03 -0.6457270 0.6610362 1
UBP4-MP5 1.233132e-02 -0.6410503 0.6657130 1
PL1-MP8 -2.571232e-03 -0.6559529 0.6508104 1
PL2-MP8 -4.812171e-03 -0.6581938 0.6485695 1
PL3-MP8 -8.325864e-03 -0.6617075 0.6450558 1
UBP1-MP8 -8.390633e-03 -0.6617723 0.6449910 1
UBP2-MP8 -5.220442e-03 -0.6586021 0.6481612 1
UBP3-MP8 -5.561452e-04 -0.6539378 0.6528255 1
UBP4-MP8 4.120582e-03 -0.6492610 0.6575022 1
PL2-PL1 -2.240938e-03 -0.6556226 0.6511407 1
PL3-PL1 -5.754632e-03 -0.6591363 0.6476270 1
UBP1-PL1 -5.819401e-03 -0.6592010 0.6475622 1
UBP2-PL1 -2.649209e-03 -0.6560308 0.6507324 1
UBP3-PL1 2.015087e-03 -0.6513665 0.6553967 1
UBP4-PL1 6.691815e-03 -0.6466898 0.6600734 1
PL3-PL2 -3.513693e-03 -0.6568953 0.6498679 1
UBP1-PL2 -3.578463e-03 -0.6569601 0.6498032 1
UBP2-PL2 -4.082710e-04 -0.6537899 0.6529734 1
UBP3-PL2 4.256025e-03 -0.6491256 0.6576377 1
UBP4-PL2 8.932753e-03 -0.6444489 0.6623144 1
UBP1-PL3 -6.476905e-05 -0.6534464 0.6533169 1
UBP2-PL3 3.105422e-03 -0.6502762 0.6564871 1
UBP3-PL3 7.769719e-03 -0.6456119 0.6611514 1
UBP4-PL3 1.244645e-02 -0.6409352 0.6658281 1
UBP2-UBP1 3.170192e-03 -0.6502114 0.6565518 1
UBP3-UBP1 7.834488e-03 -0.6455471 0.6612161 1
UBP4-UBP1 1.251122e-02 -0.6408704 0.6658928 1
UBP3-UBP2 4.664296e-03 -0.6487173 0.6580459 1
UBP4-UBP2 9.341024e-03 -0.6440406 0.6627227 1
UBP4-UBP3 4.676728e-03 -0.6487049 0.6580584 1
# Extract unique pond names from the dataframe
short_names <- unique(Mod_YF85$pond_name)
# Initialize an empty dataframe to store monthly mean results
monthly_means <- data.frame(Pond = character(), Month = character(), AvgTemp = numeric(), Year = numeric(), stringsAsFactors = FALSE)
# Convert 'time' column to Date class
Mod_YF85$time <- as.Date(Mod_YF85$time, format="%Y-%m-%d")
# Extract year and month from 'time' column
Mod_YF85$Year <- format(Mod_YF85$time, "%Y")
Mod_YF85$Month <- format(Mod_YF85$time, "%B")
# Loop through each unique pond name
for (short_name in short_names) {
# Filter data for the current pond
pond_data <- subset(Mod_YF85, pond_name == short_name)
# Loop through each month
for (month in month.name) {
# Filter data for the current month
filtered_data <- subset(pond_data, Month == month)
# Calculate the mean of 'med' for the current month
avg_med <- mean(filtered_data$med, na.rm = TRUE)
# Append results to monthly_means dataframe
monthly_means <- rbind(monthly_means, data.frame(Pond = short_name, Month = month, AvgTemp = avg_med, Year = 2020))
}
}
# Print the monthly means dataframe
print(monthly_means)
# Plotting the data
ggplot(monthly_means, aes(x = Month, y = AvgTemp, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
scale_x_discrete(limits = month.name) +
labs(title = "Monthly Average Temperatures for Each Pond",
x = "Month",
y = "Average Temperature",
color = "Pond") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Define the correct order for months
month_order <- c("January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December")
# Recalculate monthly averages with correct month ordering
calculate_monthly_averages <- function(data, decade_start, decade_end) {
# Initialize empty dataframe to store monthly averages
monthly_avg_table <- data.frame(Pond = character(), Month = factor(), AvgTemp = numeric(), Year = integer(), stringsAsFactors = FALSE)
# Loop through each year in the decade
for (yr in decade_start:decade_end) {
# Filter data for the current year
yearly_data <- data %>%
filter(year(time) == yr)
# Loop through each pond
for (pond_name in unique(data$pond_name)) {
pond_data <- yearly_data %>%
filter(pond_name == pond_name)
# Loop through each month
for (month in month_order) {
# Filter data for the current month
monthly_data <- pond_data %>%
filter(format(time, "%B") == month)
# Calculate the mean of 'med' for the current month
avg_temp <- mean(monthly_data$med, na.rm = TRUE)
# Append results to monthly_avg_table dataframe
monthly_avg_table <- rbind(monthly_avg_table, data.frame(Pond = pond_name, Month = factor(month, levels = month_order), AvgTemp = avg_temp, Year = yr))
}
}
}
return(monthly_avg_table)
}
# Calculate monthly averages for each decade
monthly_avgX50_table_2030s <- calculate_monthly_averages(Mod_YF85, 2030, 2039)
monthly_avgX50_table_2060s <- calculate_monthly_averages(Mod_YF85, 2060, 2069)
monthly_avgX50_table_2090s <- calculate_monthly_averages(Mod_YF85, 2090, 2099)
# Add a Decade column to each dataframe
monthly_means <- monthly_means %>%
mutate(Decade = "2020s")
monthly_avgX50_table_2030s <- monthly_avgX50_table_2030s %>%
mutate(Decade = "2030s")
monthly_avgX50_table_2060s <- monthly_avgX50_table_2060s %>%
mutate(Decade = "2060s")
monthly_avgX50_table_2090s <- monthly_avgX50_table_2090s %>%
mutate(Decade = "2090s")
# Combine all monthly average tables into one dataframe
combined_monthly_avg <- bind_rows(
monthly_means,
monthly_avgX50_table_2030s,
monthly_avgX50_table_2060s,
monthly_avgX50_table_2090s
)
# Print the combined dataframe
print(combined_monthly_avg)
# Define the reference year
reference_year <- 2020
# Check if the column exists and is correctly populated
summary(combined_monthly_avg)
Pond Month AvgTemp Year Decade
Length:4092 Length:4092 Min. :-0.23 Min. :2020 Length:4092
Class :character Class :character 1st Qu.: 7.42 1st Qu.:2036 Class :character
Mode :character Mode :character Median :11.18 Median :2064 Mode :character
Mean :11.18 Mean :2063
3rd Qu.:14.94 3rd Qu.:2092
Max. :21.28 Max. :2099
# Calculate average temperature for each month in the reference year
ref_year_avg <- combined_monthly_avg %>%
filter(Year == reference_year) %>%
group_by(Month) %>%
summarise(avg_temp_ref_year = mean(AvgTemp, na.rm = TRUE))
# Join reference year averages with the main dataset
combined_monthly_avg_with_diff <- combined_monthly_avg %>%
left_join(ref_year_avg, by = "Month") %>%
mutate(diff_AvgTemp = AvgTemp - avg_temp_ref_year) %>%
select(-avg_temp_ref_year) # Remove the reference year average column
# Display the structure to verify
str(combined_monthly_avg_with_diff)
'data.frame': 4092 obs. of 6 variables:
$ Pond : chr "MP1" "MP1" "MP1" "MP1" ...
$ Month : chr "January" "February" "March" "April" ...
$ AvgTemp : num 6.53 5.35 5.41 6.97 9.94 ...
$ Year : num 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
$ Decade : chr "2020s" "2020s" "2020s" "2020s" ...
$ diff_AvgTemp: num 0.487 0.49 0.482 0.487 0.504 ...
# Filter out data for the reference year
DeltaMonthly_YF85 <- combined_monthly_avg_with_diff %>%
filter(Year != reference_year)
# Display the filtered data
print(DeltaMonthly_YF85)
NA
library(dplyr)
DeltaMonthly_YF85
DiffMonthly_YF85 <- DeltaMonthly_YF85 %>%
select(Pond, Month, Year, diff_AvgTemp)
# Create a new column based on year ranges
DiffMonthly_YF85$Decade <- ifelse(DiffMonthly_YF85$Year >= 2030 & DiffMonthly_YF85$Year < 2040, "2030s",
ifelse(DiffMonthly_YF85$Year >= 2060 & DiffMonthly_YF85$Year < 2070, "2060s",
ifelse(DiffMonthly_YF85$Year >= 2090 & DiffMonthly_YF85$Year < 2100, "2090s",
"Other")))
# Print the updated data frame
print(DiffMonthly_YF85)
# Convert Month to Date type
DiffMonthly_YF85 <- DiffMonthly_YF85 %>%
mutate(Month = as.Date(paste0(Month, " 1 ", Year), format = "%B %d %Y"))
print(DiffMonthly_YF85)
# Group by Pond and Month, then calculate mean and standard deviation of diff_AvgTemp
summary_stats <- DiffMonthly_YF85 %>%
dplyr::mutate(Month = format(Month, "%m")) %>% # Extracts YYYY-MM format
dplyr::group_by(Pond, Month, Decade) %>%
dplyr::summarise(
mean_diff_AvgTemp = mean(diff_AvgTemp, na.rm = TRUE),
sd_diff_AvgTemp = sd(diff_AvgTemp, na.rm = TRUE)
)
`summarise()` has grouped output by 'Pond', 'Month'. You can override using the `.groups` argument.
# View the summary statistics
print(summary_stats)
# Group by just Month, then calculate mean and standard deviation of diff_AvgTemp
summary_stats2 <- DiffMonthly_YF85 %>%
dplyr::mutate(Month = format(Month, "%m")) %>% # Extracts YYYY-MM format
dplyr::group_by(Month, Decade) %>%
dplyr::summarise(
mean_diff_AvgTemp = mean(diff_AvgTemp, na.rm = TRUE),
sd_diff_AvgTemp = sd(diff_AvgTemp, na.rm = TRUE)
)
`summarise()` has grouped output by 'Month'. You can override using the `.groups` argument.
# View the summary statistics
print(summary_stats2)
Plotting the monthly temperature change graphs
2030s
summary_stats # grouped by month and pond
summary_stats2 # grouped by just month
# 2. Monthly averages across ponds
MonthlyStats_2030 <- summary_stats2 %>%
filter(Decade == "2030s")
ggplot(data = MonthlyStats_2030, aes(x = Month)) +
geom_point(aes(y = mean_diff_AvgTemp)) +
theme_classic()
# Convert Month from character to numeric to ensure proper ordering in ggplot
MonthlyStats_2030$Month <- as.numeric(MonthlyStats_2030$Month)
# Plot using ggplot
MonthlyStats_2030Plot <- ggplot(data = MonthlyStats_2030, aes(x = Month)) +
geom_line(aes(y = mean_diff_AvgTemp)) + # Point plot with mean values
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Horizontal dashed line at y = 0
geom_ribbon(aes(ymin = mean_diff_AvgTemp - sd_diff_AvgTemp, ymax = mean_diff_AvgTemp + sd_diff_AvgTemp), fill = "blue", alpha = 0.3) + # Ribbon for positive and negative error based on sd values
scale_x_continuous(breaks = 1:12, labels = month.abb) + # Adjust x-axis labels to show month abbreviations
ylim(-4,6)+
theme_classic() +
labs(x = "", y = "Delta Temp (C)") +
theme(axis.text.x = element_text(angle = 65, hjust = 1)) +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14))
MonthlyStats_2030Plot
2060s
# 2. Monthly averages across ponds
MonthlyStats_2060 <- summary_stats2 %>%
filter(Decade == "2060s")
ggplot(data = MonthlyStats_2060, aes(x = Month)) +
geom_point(aes(y = mean_diff_AvgTemp)) +
theme_classic()
# Convert Month from character to numeric to ensure proper ordering in ggplot
MonthlyStats_2060$Month <- as.numeric(MonthlyStats_2060$Month)
# Plot using ggplot
MonthlyStats_2060Plot <- ggplot(data = MonthlyStats_2060, aes(x = Month)) +
geom_line(aes(y = mean_diff_AvgTemp)) + # Point plot with mean values
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Horizontal dashed line at y = 0
geom_ribbon(aes(ymin = mean_diff_AvgTemp - sd_diff_AvgTemp, ymax = mean_diff_AvgTemp + sd_diff_AvgTemp), fill = "blue", alpha = 0.3) + # Ribbon for positive and negative error based on sd values
scale_x_continuous(breaks = 1:12, labels = month.abb) + # Adjust x-axis labels to show month abbreviations
ylim(-4,6)+
theme_classic() +
labs(x = "", y = "") +
theme(axis.text.x = element_text(angle = 65, hjust = 1)) +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14))+
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank())
MonthlyStats_2060Plot
2090s
# 2. Monthly averages across ponds
MonthlyStats_2090 <- summary_stats2 %>%
filter(Decade == "2090s")
ggplot(data = MonthlyStats_2090, aes(x = Month)) +
geom_point(aes(y = mean_diff_AvgTemp)) +
theme_classic()
# Convert Month from character to numeric to ensure proper ordering in ggplot
MonthlyStats_2090$Month <- as.numeric(MonthlyStats_2090$Month)
# Plot using ggplot
MonthlyStats_2090Plot <- ggplot(data = MonthlyStats_2090, aes(x = Month)) +
geom_line(aes(y = mean_diff_AvgTemp)) + # Point plot with mean values
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Horizontal dashed line at y = 0
geom_ribbon(aes(ymin = mean_diff_AvgTemp - sd_diff_AvgTemp, ymax = mean_diff_AvgTemp + sd_diff_AvgTemp), fill = "blue", alpha = 0.3) + # Ribbon for positive and negative error based on sd values
scale_x_continuous(breaks = 1:12, labels = month.abb) + # Adjust x-axis labels to show month abbreviations
ylim(-4,6)+
theme_classic() +
theme(axis.text.x = element_text(angle = 65, hjust = 1)) +
labs(x = "", y = "") +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14))+
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank())
MonthlyStats_2090Plot
Combine the plots together here for now
# Combine plots into a panel plot
DeltaMonth_YF85_plot <- ggarrange(MonthlyStats_2030Plot, MonthlyStats_2060Plot, MonthlyStats_2090Plot, ncol = 3)
# Save the combined plot
ggsave("Corr_DeltaMonthYF85.png", plot = DeltaMonth_YF85_plot, width = 12, height = 3)
Create .csv files of the data for these graphs
write.csv(MonthlyStats_2030, file = "Corr_MonthlyStatsYF85_2030.csv")
write.csv(MonthlyStats_2060, file = "Corr_MonthlyStatsYF85_2060.csv")
write.csv(MonthlyStats_2090, file = "Corr_MonthlyStatsYF85_2090.csv")